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Abstract 


A numerical study of the steady, axi sym- 
metric flow in a heated, rotating spherical 
shell is conducted to model the Atmospheric 
General Circulation Experiment (AGCE) proposed 
to run aboard a later shuttle mission. The AGCE 
will consist of concentric rotating spheres con- 
fining a dielectric fluid. By imposing a di - 
electric field across the fluid a radial body 
force will be created. The numerical solution 
technique is based on the incompressible Navier- 
Stokes equations. In the method a pseudospec- 
tral technique is used in the latitudinal dir- 
ection, and a second-order accurate finite dif- 
ference scheme discretizes time and radial der- 
ivatives. This paper discusses the development 
and performance of this numerical scheme for the 
AGCE which has been modelled in the past only by 
pure FD formulations. In addition, previous 
models have not investigated the effect of using 
a dielectric force to simulate terrestrial 
gravity. The effect of this dielectric force on 
the flow field is investigated as well as a 
parameter study of varying rotation rates and 
boundary temperatures. Among the effects noted 
are the production of larger velocities and 
enhanced reversals of radial temperature 
gradients for a body force generated by the 
electric field. 
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where g is the acceleration due to gravity, a is 
the thermal expansivity, AT is the horizontal 
temperature difference, d is the annulus height, 
(r 0 - r i ) is the radial thickness, !J the 
rotation rate, and u the kinematic viscosity. 
The knee-shaped curve in Figure 1. referred to 
as the critical curve, separates the axis- 
ymmetric flow regime from the wave-like flow 
(nonaxi symmetric regime). 

The flow in the nonaxisymmetric regime is 
an analogue of the synoptic scale motion of a 
planetary atmosphere. The annulus experiments 
have therefore guided studies concerning the 
predictability or periodicity of the atmos- 
phere. Scaling of the governing equations based 
on large scale features of the atmosphere and 
laboratory experiments leads to the same set of 
approximate equations for both systems. Thus, 
the experiments are models of the synoptic-scale 
atmospheric flow. 


I. Background 

Since as early as the nineteenth century 
there have been experimental efforts aimed at 
better understanding the earth's atmosphere. 
Until the 1950's, however, these efforts were 
qualitative in nature. It was found that exper- 
iments involving a fluid confined between two 
concentric, rotating cylinders could be syste- 
matically studied. Experiments performed by 
Hide 1 in 1953, using a cylindrical annulus, 
resulted in the discovery of a meandering jet 
which greatly resembled the earth's jet stream. 

In 1965, Fowl is and Hide 2 quantitatively 
produced a flow regime diagram by careful 
measurements taken from the annulus 
experiments. Figure 1. is the regime diagram 

adapted from Fowlis and Hide. The ordinate is a 
thermal Rossby number, and the abscissa a Taylor 

number, both of which are given in equations (1) 

and (2): 


Since analytical solutions are not obtain- 
able due to the strong nonlinearity of the gov- 
erning equations, numerical modelling has served 
to enhance the theory and provide a quantitative 
comparison with experiment. Although models are 
becoming increasingly more complex, substantial 
advances in atmospheric studies have come from 
simple models which have isolated essential 
mechanisms concerning the earth's general circu- 
lation. It was observed by Lorenz 2 * * * that a 
steady axisymmetric solution to the governing 
equations exists everywhere in parameter space. 
The reason for the nonaxisymmetric regimes is 
the instability of the axisymmetric flow to 
infinitessimal (nonaxisymmetric) pertubations. 
The first step, therefore, in determining 
whether a given flow is stable is to calculate 
the steady, axisymmetric states. If a given 
axisymmetric solution is unstable to perturba- 
tions of any wavelength, then it might not be 
observed in a physical situation. Axisymmetric 
states for the thermally driven motion of a 
fluid contained in a rotating cylindrical 
annulus were numerically investigated by 
Will iams4«5. Axisymmetric states were 
thoroughly studied, but a variety of points in 
parameter space were not considered. Miller and 
Gall 6 used a linear model to predict a tran- 
sition curve very close to that observed in the 
laboratory. 
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An obvious disadvantage of cylindrical 
models of the earth's atmosphere Is the 
geometrical difference. However, experimental 
studies done in a spherical geometry on earth 
cannot generate a radial body force strong 
enough to compete with terrestrial gravity. 
Such an experiment needs the near zero gravity 
of Spacelab to generate a sufficiently dominate 
radial body force. 

The numerical model in this study will 
serve as an aid in constructing a regime diagram 
for a spherical geometry, analogous to the 
cylindrical annulus case (Figure 1.) The 
spherical regime diagram is also most essential 
for the AGCE , since knowing the location of the 
axisymmetric and nonaxisymmetrlc flow regimes is 
necessary for the design of the AGCE Instrument. 

In addition to the need for a spherical 
regime diagram for AGCE, it Is desirable to know 
the effect of the radial variation of the di- 
electric gravitational field. The radial force 
created by imposing an electric field on the di- 
electric fluid is Inversely proportional to the 
fifth power of the radius’. On earth, the 
gravitational force is inversely proportional to 
the square of the radius. The effect of varying 
this power is needed to determine how capable 
the dielectric force is in simulating ter- 
restrial gravity. 


II. Physical Model 

A schematic of the experimental set-up is 
shown in Figure 2. A dielectric fluid will be 
confined in a portion of the gap between the 
spheres by placing solid walls in the spherical 
domain. The annulus is subdivided to provide 
room for instrumentation. Temperature profiles 
on the spheres will simulate the equator to pole 
distribution and establish the large-scale ver- 
tical stability (warmer outer sphere) of the 
atmosphere. The gravitational body force on the 
atmosphere is simulated by application of an 
electric field across the gap which generates a 
radial body force in the dielectric fluid. In 
order to be radial, this force requires the zero 
gravity environment of space, since the earth's 
gravitational force would greatly distort the 
simulated gravitational force field. 


and lower temperature extremes on the 
spherical walls. The nondimensional (unprimed) 
variables vorticity ( 5 ), stream function (v), 
and swirl (<a) are defined in equations (31 
through (7). 
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equations (where unprimed variables are here 
nondimensional) are: 
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III. The Numerical Model 
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The right-handed coordinate system used in 
the model is rotating with the apparatus at the 
constant angular velocity ft. Figure 3. depicts 
this system with r representing the radius, 9 
the co-latitude, and ♦ the azimuthal angle. The 
governing equations are solved in a vorticity/ 
stream function formulation for a Boussinesq 
fluid. The solid body rotation rate, fl, non- 
dimensionalizes time. The characteristic radial 
velocity is taken as the solid body rotation 
rate times the gap width— all other length 
scales are nondimensionalized by the outer 
sphere radius. A nondimensional temperature is 
obtained from the difference between the upper 
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Boundary conditions appropriate for the AGCE 
were Incorporated, i.e., the no-slip conditions 
and conducting walls. Details of the numerical 
model may be found in previous work®. A brief 
summary of the scheme is now given. 

As stated earlier the numerical model 
utilizes a mixed pseudospectral /finite differ- 
ence (PS/FD) scheme. Due to the presence of the 
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latitudinal side walls (see figure 2.), 
Tschebyshev polynomials are used as the 
spectral expansion functions. In addition 
Tschebyshev collocation methods insure an ex- 
ponential rate of convergence as opposed to 
the algebraic rates associated with FD 
methods 9 . It was sufficient to use the 
pseudospectral formulation In only one direction 
(latitude) since the computational grid required 
1/7 to 1/3 the number of grid points as the pure 
FO model . 

The numerical scheme treats the rotation 
and advection terms explicity while radial 
diffusion terms are implicity solved. 
Pseudospectral derivatives are computed by 
setting up appropriate matrices for Tschebyshev 
collocation outside the time iteration process. 
Latitudinal derivatives are thus computed via a 
simple matrix multiplication at each time step. 

The Poisson equation for stream function 
was solved in two ways. Initially a pseudo-time 
derivative was introduced so that at steady 
state the Poisson equation would be satisfied. 
During this internal iteration latitudinal der- 
ivatives were not updated, which introduced a 
stabilizing factor in the scheme. This effec- 
tive "pseudo-time lag" in the Poisson iteration 
allowed the model to utilize a time step compar- 
able to an implicit formulation in the overall 
iteration scheme. Alternatively, the matrix 
corresponding to the Laplacian operator was in- 
verted by a conventional LU decomposition and 
back solution. Though a fewer number of steps 
were required in the latter formulation a much 
more stringent time Interval restriction was 
found necessary. 


Numerical Model Application: 
Experimental Parameters 


The main results of an extensive parameter 
study which varied boundary temperatures and ro- 
tation rates for several gravitational force 
distributions will now follow. The latitudinal 
side walls were set at ir/4 and 3n/4. The inner 
and outer sphere radii were set at 5 and 6 cm., 
respectively. The term governing the force 
distribution is the buoyancy term of the 
vorticity equation (Eq. 8) as given below: 
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Cases I through IV, defined in Table 
1., allow the study of the effect on the flow 
field if Gr and Ek are allowed to vary by four 
orders of magnitude while maintaining a uniform 
gravitational field (n = 0). Each of these 
quantities have physical significance. The 

Grashof number represents the ratio of buoyancy 
to viscous forces. It appears in the buoyancy 
term of the vorticity equation numbered 8. The 
Ekman number measures the ratio of viscous to 
Coriolis forces in the flow field. This 

parameter is found in the vorticity and energy 
equations. If Gr * 0, it couples each of the 

other equations through the vorticity equation. 
If Gr = 0, the energy equation becomes 

uncoupled. 

Figures 4. through 9. are contour plots of 
w, u, v, T, and 5 for all four cases. Most 
noted in the contour plots is the increasing 
departure from symmetry. In addition, increased 
meridional transport for Increasing Gr is 

evident in plots of stream function (Figure 4.) 
The radial and zonal jets seen in the contours 
of w (Figure 5.) and u (Figure 6.), 
respectively, are seen to move northward as 
buoyancy forces Increase (higher Gr). The 

latitudinal jets associated with v (Figure 7.) 
become more concentrated towards the spherical 
wall. Strong thermal boundary layers develop as 
is evident in temperature contour plots of Case 
IV, relative to the previous cases, which for 
Case I is essentially a conduction profile 
(Figure 8.). Vorticity maxima are seen in 
Figure 9. to shift from the inner radial wall at 
mid-latitudes for Case I to the wall at */4 and 
approximately 5.6 cm for Case IV. These effects 
are all associated with an Increasing buoyancy 
force which appears to be the dominant factor 
relative to increasing rotation by the same 
order of magnitude. 

Listed in Table 1. is the product GrEk2 
calculated for the four cases. This product re- 
presents the ratio of buoyancy to Coriolis for- 
ces. GrEk? is the smallest in Case IV 

which has the largest flow field values. The 
enhancement of advection for decreasing values 
of GrEk 2 can be explained as follows. With a 
radial gravitational field, warm fluid will tend 
to rise toward the outer sphere. For large 
GrEk 2 numbers, buoyancy forces more 
effectively keep warm fluid near the outer 
sphere. If buoyancy forces are weakened, 
meridional flow will increase to enforce the 
stable stratification of fluid in the interior. 


An interesting study regarding heat trans- 
In the cases to be discussed, the exponent, n, fer rates for the four cases was made by carry - 

and gravitational acceleration, g, were varied ing out a heat balance on the domain from */4 to 

under the restriction that the mean gravita- x/2 and from r-j to r 0 . Ratios of the amount 

tional acceleration be maintained at 980 cm/s 2 of heat transfer occurring at a given boundary 

This was done by calculating a g re f as given over that occurring in pure conduction are 

in (11): 



3 



listed in Table 1. The temperature isotherms of 
Figure 8c. (Case III) and 8d. (Case IV) indicate 
negative radial thermal gradients around 6 » 
5*/12 near the lower spherical wall. This is 
due to the increasing meridional motion for 
these two cases relative to Cases I and II which 
carry warm fluid northward as cold fluid sinks 
towards the equator. This motion acts to- 
decrease the radial thermal gradients near the 
spherical walls and, for cases of strong 
meridional motion, to even reverse these 
gradients. The result is a decrease in the 
amount of heat transfer occurring at the 
spherical walls as evidenced in Table 1. 

It should be noted that at steady state no 
net torque is exerted on the fluid. Thus the 
sum of the viscous stresses exerted by the 
boundaries on the fluid should vanish. This 
means that the distribution of easterly and 
westerly flows at the boundaries should be about 
equal throughout the fluid volume. Such a trend 
is seen in the contour plots of azimuthal 
velocity (Figure 7). However, the distribution 
of easterlies and westerlies in the atmosphere 
is different, which is expected since there Is 
no outer boundary as is the case for the AGCE. 
Time-averaged atmospheric data yield the polar 
easterlies, the mid-latitude westerlies and the 
tropical easterly trade winds. In the AGCE 
model, contour plots of azimuthal velocity 
reveal a region of positive zonal flow 
(westerlies) in the upper half of the domain, 
while the lower portion, down to the equator 
shows negative zonal velocity (easterlies). 
These flows may be studied using conservation 
principles. If a particle at the equator, 
having initially no zonal (easterly or westerly) 
component in its motion, goes toward the polar 
axis while conserving its angular momentum, it 
is flung eastward. Conversely, a particle going 
toward the equator from one of the poles would 
have a deficit of west-to-east momentum with 
respect to the earth and would have to curve 
toward the west. The resulting zonal linear 
speeds with respect to the earth's surface can 
be computed by applying the principle of 
conservation of angular momentum for a unit 
mass. However, the observed wind speeds in the 
atmosphere are substantially less than those 
calculated on the basis of angular momentum 
conservation. It is believed that in the 
earth's atmosphere the meridional flow breaks 
down into smaller circulations or cells. 10 
Throughout the atmosphere the meridional 
components are on the average less than 
one-tenth of the zonal components. It was found 
in this numerical study that meridional flow 
values were on the average one-third of the 
zonal flow. The break-up of the meridional flow 
into two cells was noted to occur under specific 
conditions. At approximately the same location 
where temperature contours Indicated a gradient 
reversal, the stream function contours were 
found to have two smaller cells in addition to 
the large cell representing the main meridional 
circulation. The reversal of the latitudinal 
gradient of temperature is certain to be a 
factor in causing these recirculations. It was 
noted that these recirculations were stronger 
for values of n equal to 0 and 2. Higher values 
of n were also found to increase flow field 
values. It is possible that these higher values 


are due to the weakening of these recirculating 
cells in the meridional flow for the cases of n 
= 5. Stream function contours show the appear- 
ance of these recirculations which are seen to 
be more prevalent for n » 0 as opposed to n « 2 
and even less pronounced when n » 5. Stream 
function contours where rotation was increased 
from 1 rad/s to 3 rad/s showed the appearance of 
these smaller cells for increased rotation, 
which also Inhibits meridional flow. 


VI. Discussion 

Experimentation with the dielectric force 
(n * 5) should be observed with the knowledge 
that differences do exist between the resulting 
flow field and a flow calculated with a force 
proportional to l/r‘ (the earth's gravita- 
tional force). Flow field values will tend to 
be higher for the dielectric force. The 
meridional flow's breakup into two smaller 
cells, which has been observed in the atmosphere 
and shown numerically for the cases of n « 0 and 
n « 2, will be less pronounced for flows 
generated with a dielectric force. To compen- 
sate for these differences, rotation rates 
should be at least greater than I rad/s since 
meridional flow is Inhibited by rotation. In 
addition, boundary temperature differences 
should not be too extreme so that Internal 
temperature gradients will be weak enough to 
allow gradient reversals and subsequent recir- 
culations in the meridional flow to occur. 
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Table 1. Heat Transfer Ratios for Increasing Grashof and Ekman Numbers. 


Case 

Gr 



Ek 


GrEk 3 

q */4 

q */2 

q 5cm 

q 6cm 

q (it/4)C* 

q {f/2)C 

q 5C 

q 6C 

I 

2.84 x 

10 3 

2.78 

X 

lO" 1 

219.2 

1.0 

0.0 

1.0 

1.0 

II 

2.84 x 

10 4 

2.78 

X 

io - 4 

21.92 

2.0 

0.0 

0.861 

0.913 

III 

2.84 x 

10 s 

2.78 

X 

10* 3 

2.192 

8.8 

0.0 

0.566 

0.890 

IV 

2.84 x 

10 6 

2.78 

X 

10- 4 

0.2192 

15.8 

0.0 

0.287 

0.898 


*q c = heat transfer from conduction solution 



Fig. 1. Cylindrical Annulus Regime Diagram. 
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Fig. 4. :p Contour Plots for Cases I to IV. 



Fig. 5. w Contour Plots for Cases I to IV. 
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Fig - 7. V Contour Plots for Cases I to IV. 




Fig. 8. T Contour Plots for Cases I to IV. 



Fig. 9. 5 Contour Plots for Cases I to IV. 

a ’ ^ max 3 •00041s"* at 1.178 rad, 5 cm. 
b ‘ ^ max = -OOSSs" 1 at 1.178 rad, 5 cm. 
c ' S max = ,324s-1 at -7854 rad, 5.5 
d - 5 max 3 i- 7195 * 1 at *7854 rad, 5.617 cm. 
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